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Computer Oriented Programs on Multi-Dimensional 

Numerical Integration 

D.A.GISMALLA & Dr. F. M. DAWALBAIT 


Abstract — Computer Oriented Programs are SOFTWARE 
BULIT-IN FUNCTIONS OR WRITTEN SOFTWARE in a 
different COMPUTER LANGUAGES. In the past only 
NUMERICAL SOFTWARE ROUTINES available are NAG 
ROUTINES written in FORTRAN and ALGOLW which both 
are POCESSING LANGUAGES and then later C++ languages 
.The FORTRAN LANGUAGE can compute to a very high 
accuracy but require and need MAIN-FRAMES 
COMPUTERS . Instead people now a days are inclined to use 
PERSONEL COMPUTERS for the money cost is few in such 
away any person can get one . 

In this paper We FIRST list some of the 
FUNCTIONS BULIT IN MATHEMATICA OR MATLAB 
programming languages that deals and compute numerically or 
( analytically in closed forms for 1-dimentional ) and higher 
dimensional integrals. 

SECOND, We will design a MATLAB SOTWARE 
PROGRAMS for 1-dimentional problems or higher 
multi-dimensional. For 1-dimentional Integrals METHODS We 
use are Trapezoidal, Simpson’s and Gaussian Quadrature 
Rules. For 2-dimentional Integrals We apply Radon’s approach 
and Gismlla’s approach with regions are the SQUARE 
REGION and RECTANGLE REGION WITH WEGIHT 
FUNCTION THE SQUARE ROOT of X, respectively .THIRD, 
for 3-dimentional or higher Integrals We apply Levin’s 
Transform with our expectation to evaluate Integrals to a high 
accuracy ((as We have dealt with before in[ 4] but using 
FORTRAN languages on main frame COMPUTERS )). 

Further, the reader can be acquainted with some symbolic 
languages and distinguishes them from the usual processing 
languages. Furthermore, the reader will know that LEVIN’S 
TRANSFORM can be used as one of the best method to sum a 
large class of series to a high accuracy. 

Index Terms — Symbolic Languages and Built-in Functions 
.Simpson’s, Trapezoidal, Gaussian, Radon’s, Gismalla’s and 
Levin’s Rules. Full Machines Accuracy and Main Frame 
Computers. 


I. BACKGROUND AND REVIEWS 

The problem of numerical integration is one of the oldest 
problem in mathematics that gives a beautiful and an insight 
light to a wide range of theoretical and computational 
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Methods in Mathematics and Numerical Analysis. This can be 
seen when literature retrieval for theory and algorithms is 
seek, especially for one-dimensional Quadrature Formulae 
such as Lobbotta, Newton’s, Trapezoidal, Simpson and 
Gaussian Quadrature Rules or for two-dimensional Cubature 
Rules as in the references [6] ,[8] ,[12],[17],[18], [19] and 
[ 20 ]. 

The theory of Cubature is burst out when the Product 
Rules generated from the one-dimensional Rules such as 
Simpson Rules having the disadvantages of accumulating 
rounding errors and do not attain the required minimum 
number of points for functions evaluations in the formula. 
Stroud [17] has developing many cubature formulae on 
symmetric regions such as square and the cubic. He creates 
by theoretical approaches to formulate a matrix system of 
equations and then solve this system by choosing and 
selecting a special paten of points to the required cubature 
formula. The problem of selecting such paten of points in 
such away to solve the matrix system is difficult unless a 
particular suitable paten is chosen . This can be seen as a 
disadvantage. However, if these suitable paten of points can 
be chosen in advance in such away to solve the system of 
equations , then many cubature formulas can be obtained. 

Cohen & Gismalla [6] has selected a certain paten 
of points on the square and the cubic to construct some types 
of cubature formulas known as symmetric cubature formulas 
as in [6 ]. Selecting in advance the suitable paten of points can 
construct good formulas , not for symmetric regions only but 
on any region with a weight function for that region as in 
Gismalla [7]. 

The Germanium ( or the Russian ) Radon [14] 
investigates and develop a large information theories on 
polynomials , orthogonal polynomials and zeros polynomials 
to establish his fifth degree cubature formula on the square 
region. If someone reads the translation for this work, he will 
certainly acknowledge the valuable efforts done to connect 
these theories in such away to establish Radon Rule. One of 
the neatness and beauty of Radon approach is that Gismalla [ 
7] has summarized it in steps as an algorithm in such away to 
construct cubature formulas not on a symmetric region as a 
square but on any region with its weight function . This can 
be seen in Gismalla[7 ] for deriving the same fifth degree 
formula on rectangle region 

£> = ;i] < x < 1 end - 1 < x < 1} with weight function 
V x . 

The theories of polynomials as investigated in Radon 
Rule as algorithms to construct cubature formulas was carried 
out further in Moller [23] and Shimds [24 ]. Moller seeks 
formulas that attain the minimum odd number of points while 
Shimd’s work[24 ] to attain the minimum even number of 
points . Gismalla[8] apply Schimd’s approach to construct a 
forth degree formula and two sixth degree formula on the 
rectangle region £1 = [0 < x < 1 and — 1 < x < 1} with 
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weight function \fx . Both these two approaches can be 
computerized to generate many cubature formulas but the cost 
of complexity is too great and high . Even if the complexity is 
revealed one can expect to obtain formulas with some points 
outside the region in consideration. 

Nevertheless, Stroud [16] has written a great deal of 
MATLAB SOFTWARE for many regions such as the circle, 
square, cubic and triangle. The useful of formulas on a 
triangle can be dealt with for problems on finite element 
methods especially when it is symmetric. NAG 
ROUTINES in FORTRAN languages are available but with 
the restriction that most of them on main frames computers. 
QUADPACK for numerical Integration can found on 
FORTRAN or C++ languages or elsewhere. 


II. THE BULIT-IN MATHEMATICA INTEGRATE 

FUNCTIONS 

Here , We discuss some of the built_in symbolic Mathematica 
functions for integration . These functions are Integrate , 
NIntegrate and Manipulate with Integrate. 

A. Integrate 

The Integrate function is used to evaluate the 

jn 1 

indefinite integral I dx 

in Mathematica Mode environment . Typing just the 
command Integrate in the following line below and then press 
the two keys SHIFT+ENTER simultaneously to get the result 
for required indefinite integral as in Example 2.1. The reader 
must know a little knowledge about how to execute a 
Mathematica function command . Once the command 
Integrate[l/(x A 3+l),x] is typed and the two keys 
SHIFT+ENTER are pressed simultaneously , the 
Mathematica will execute it as input automatically with 
in[l]:= to indicate that this the first input as 

in[l]:= Integrate[l/(x A 3+l),x] 

The output result will be preceded automatically with out[l] := 
to indicate that this the first output result as an answer 

AreTnq ■ 1 _!~ i i 

out[l]:= -—+ - Log[l+x]- - Log[l-x+x 2 ] 

Similarly , any further input or output will be preceded by its 
number in sequence automatically for second input or output 

in[2]:= or out[2]:= ,.... etc. 

Example 2.1 

The following Integrate function gives the indefinite 
integrals 

Integrate \f,x\ , gives the indefinite integral j f d x 

In[l]:= Integrate[l/(x A 3+l),x] 

t i , 

Out[l]= - - —-1- - Log[l+x]- - Log[l-x+x 2 ] 


In[4]:= Integrate [Log [l-x A 2],x] 

Out[4]= -2x - Log[-l+x] + Log[l+x] + x Log[l-x 2 ] 

The following Integrate function gives the definite 

integral 

Integrat e\f 9 {x^c min ^c max }] 4he definite 


integral f(#)dx 

Here are Integrate functions for the definite 
integral J E f(jc)dx . 

In[5] := Integrate[Sin[x] A 2, {x,a,b} ] 
Out[5]= ~(- a + b + Cos[a]Sin[a] - Cos[b]Sin[b]) 


In[6]: = Integrate [Exp [-x A 2], {x,0,Infinity} ] 


-'Tl 


Out[6]= y 


In[7]:= Integrate[l/(x A 3+l),{x,0,l}] 

Out[7]= ~ (2\ 3 n +Log[ 641) 


In[8]:= Jp' 

Out[8]= — -n:(EulerGamma+Log[4]) 


Mathematica cannot give you a formula for this definite 
integral Jp dx but instead We can get a numerical 

result 

In[9]:= Integrate [x A x],{x,0,1}] 

Out[9]= j $x x dx 


Here below N and % indicates the numerical for the current 

input % Integrate Command 
In[10]:= N[%] 

Out[10]= 0.783431 


Here is the Integrate for the multiple integral 


r™-I:™-flyldy 


O'TKJJC 

"TteJtl * " V-min 

Integratef/^lx^x^Tj^x^^Jjly •>ymm , >ymax\ 


For Multiple integral with x integration outermost: 
In[ll]:= Integrate[Sin[x y],{x,0,l },{y,0,x}] 

Out[ll]= 7 (EulerGamma-CosIntegral[l]) 

u 


In[12]:= J 0 L Shi [acy] dy dac 
Out[12]= ~ (EulerGamma-CosIntegral[l]) 

■ill 

Integrals over Regions 

This does an integral over the interior of the unit circle. 
In[13]:= Integrate [If [x A 2+y A 2<l,l,0],{x,-1,1 },{y,-l,l}] 

Out[13]= n 
Here is an equivalent form. 


In[2]:= Integrate [x A n,x] 

u+l 

Out[2]= Z— 

ia + 1 


In[3]:= 
Out[3]= - 


Integrate [ 1 /(x A 4-a A 4) ,x] 
JrcTan[ Logfj-x] Log[a+s] 
2^3 + 4i 3 4i 3 


In[14]:= Integrate[Boole[x A 2+y A 2<l],{x,-l,l },{y,-l,l}] 

Out[14]= n 

Even though an integral may be straightforward over a simple 
rectangular region, it can be significantly more complicated 
even over a circular region. 

This gives a Bessel function. 
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In [15]:= Integrate [Exp [x] 
Boole[x A 2+y A 2<l],{x,-l,l },{y,-l,l}] 
Out[15]= 2 n Bessell[l,l] 


B. NIntegrate 


NI n teg r a te I/, jc max } ] 

gives a numerical approximation to the integral 

NIntegrate\f,{xpc min pc max },{y min^y max\ i* • •] 

gives a numerical approximation to the multiple integral 


r 


This implies that the built-in NIntegrate Mathematica 
function evaluates integrals in one or multi-dimensional 
integrals numerically to a good limit of accuracy as can be 
shown here in Example 2.2 

Example2.2 Compute a numerical integral: 

In[ 16] := NIntegrate[Sin[Sin[x]], {x,0,2} ] 
Out[16]= 1.24706 


This finds a numerical approximation to the integral 



dx. 


In[ 17] := NIntegrate [Exp [-x A 3 ], {x,0,Infinity} ] 

Out[17]= 0.89298 


Compute a multi dimensional integral (with singularity at the 

origin): 

In[18]:= 

NIntegrate[-j====,{xl,0,l},{x2,0,l},{x3,0,l},{x4,0,l 

^ 3 + i + 2 _ ia 

}] 

Out[18]= 1.2391 


Here is the numerical value of the double integral 

TiTiC^ + yVy dx 


In[19]:= NIntegrate [x A 2+y A 2,{x,-1,1 },{y,-l,l}] 

Out[19]= 2.66667 


Out[20] 



In[21 ] := Manipulate[Integrate[Sin[x( 1 +a 
x)],{x,0,6}],{a,0,2}] 


Out[21]= 



In[2 2]:= Manipulate[Integrate[l/(x A n-l),x],{n,2,10,l}] 
Out[22]= 


n f - 

n 

0 

-n 

1 1 c - 

- ; u 

ArcTanLkU 1 

□-□ - Log. 

2 4 

.□1 Dx. 

1 i 

.□ - LogUL □ x_ 

4 



C. Manipulate with Integrate 

The Command Maipulate means that parameters can be 
manipulate continuously or in discrete steps. Also , it is 
possible to manipulate two parameters . The manipulation can 
be done for many results such as expansion or plotting a 
function. Here , We give an Example 2.3 for Manipulate with 
Plot and Integrate 

Example 2.3 (( See In[20] —In[22])) 


In [20]:= Manipulate[Plot[Sin[x (1+ax)], 
{x,0,6}],{a,0,2}] 


III. THE BULIT-IN MATLAB QUAD FUNCTIONS 

Here We shall describe the four Built-in Matlab QUAD 
functions to evaluate numerically an integrand over 
One-dimensional accurately to about six decimal places. 
These functions are QUAD, QUADL, DBLQUAD and 
TRIPLEQUAD as can be seen in sections 3.1, 3.2, 3.3 and 
3.4, respectively. 

A. QUAD 

QUAD Numerically evaluate integral, adaptive Simpson 
quadrature 

Q = QUAD(FUN,A,B ) 


Q = QUAD(FUN,A,B) tries to approximate the integral of 
function FUN from A to B to within an error of 1 .e-6 using 
recursive adaptive Simpson quadrature. The function Y = 
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FUN(X) should accept a vector argument X and return a 
vector result Y, the integrand evaluated at each element of X. 

Q = QUAD(FUN,A,B,TOL) uses an absolute error 
tolerance of TOL instead of the default, which is l.e-6. 
Larger values of TOL result in fewer function evaluations and 
faster computation, but less accurate results. The QUAD 
function in MATLAB 5.3 used a less reliable algorithm and a 
default tolerance of l.e-3. 

Use array operators.*, ./ and . A in in the definition of 
FUN so that it can be evaluated with a vector argument. 

jn2 i 

Example 2.4 Approximate the integral J : 3 c-r by 
using quad function in 

Matlab. Observe that the FUN in Q=QUAD(FUN , A,B) can 
be sent as an argument to the function quad using two 
approaches i.e. FUN can be specified as: 

An inline object: 

F = inline(T./(x. A 3-2*x-5)'); 

Q = quad(F,0,2); 

A function handle: 

Q = quad(@myfun,0,2); 
where myfun.m is an M-file: 
function y = myfun(x) 
y = l./(x. A 3-2*x-5); 


Type in the Command window to demonstrate an inline object 

» F = inline(T./(x. A 3-2*x-5)'); 

» Q = quad(F,0,2) 

Q = -0.4605 


Or Type in the Command window to demonstrate a function 

handle: 

» Q = quad(@myfun,0,2) 

Q= -0.4605 


B. QUADL 

The differences between the built-in function QUADL and 
recursive adaptive QUAD is that it uses high order 
Lobatto quadrature and the function QUAD may be more 
efficient with low accuracies or nonsmooth integrands . 
Further QUADL can be executed similary as QUAD for FUN 
can be specified as inline or as a function handle. 

Example 2.5 Approximate the integral J : -3i __ by 

using quadL function in Matlab 

Type in the Command window to demonstrate an inline object 

» F = inline(T./(x. A 3-2*x-5)'); 

» Q = quadl(F,0,2) 

Q= -0.4605 


Or Type in the Command window to demonstrate a function 

handle: 

» Q = quadl(@myfun,0,2) 

Q= -0.4605 


C. DBLQUAD 

DBLQUAD Numerically evaluate double integral. 
DBLQUAD(FUN,XMIN,XMAX,YMIN,YMAX) evaluates 
the double integral of FUN(X,Y) over the rectangle XMIN 
<= X <= XMAX, YMIN <= Y <= YMAX FUN(X,Y) should 
accept a vector X and a scalar Y and return a vector of values 
of the integrand . 

DBLQUAD(FUN,XMIN,XMAX,YMIN,YMAX,TOL) 
uses a tolerance TOL instead of the default, which is l.e-6 . 

DBLQUAD(FUN,XMIN,XMAX,YMIN,YMAX,TOL,@ 
QUADL) 

uses quadrature function QUADL instead of the default 
QUAD FUN can be an inline object or a function handle 

Example 2.6 Approximate the integral 

O C(y * + x * cosCy)) f£y}.djr by 

using dblquad function in Matlab 

Type in the Command window to demonstrate an inline object 

» F = inline('y*sin(x)+x*cos(y)'); 

» Q = dblquad(F, pi, 2*pi, 0, pi) 

Q = -9.8696 


Or Type in the Command window to demonstrate a function 

handle: 

» Q = dblquad(@integrnd, pi, 2*pi, 0, pi 

Q =- -9.8696 

Observe that integrnd.m is an M-file: 
function f = integrnd(x, y ) 
f = y*sin(x)+x*cos(x); 


D.TRIPLEQUAD 

TRIPLEQUAD Numerically evaluate triple integral. 

TRIPLEQUAD(FUN,XMIN,XMAX,YMIN,YMAX,ZMIN, 
ZMAX) evaluates the triple integral of FUN(X,Y,Z) over the 
three dimensional rectangular region XMIN <= X <= 
XMAX, YMIN <= Y <= YMAX, ZMIN <= Z <= ZMAX 
FUN(X,Y,Z) should accept a vector X and scalar Y and 
Z and return a vector of values of the integrand. 
TRIPLEQUAD(FUN,XMIN,XMAX,YMIN,YMAX,ZMIN, 
ZMAX,TOL) uses a tolerance TOL 
instead of the default, which is 1 .e-6 
Example 2.7 
Approximate the integral 

C[ Jn £ S-Jy *sin(xi -h z* cosiyi^dz }dy\dx 
by using triplequad function in Matlab 


Type in the Command window to demonstrate an inline object 

» F = inline('y*sin(x)+z*cos(y)'); 

» Q = triplequad( F, 0, pi, 0, 1,-1, 1) 

Q = 2.0000 


Or Type in the Command window to demonstrate a function 

handle: 

» Q = triplequad(@integrnd, 0, pi, 0, 1,-1, 1) 
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Q =- 2.0000 

Observe that integrnd.m is an M-file: 
function f = integrnd(x, y, z) 
f = y*sin(x)+z*cos(x); 


IV. MALAB SOFTWARE METHODS & ALGORITHMS 
FOR NUMERICAL INTEGRATION 


These techniques are a generalization of a small low-order 
degree formula. The main interval is piece wisely divided 
into small sub-intervals and the required rule is applied on 
each ones. If we apply Simpson's rule on each subinterval it 
is called Simpson’s composite rule as shown by Theorem 
(4.2) while Trapezoidal composite rule by the following 
Theorem (4.1) 

A. Trapezoidal Rule 


The composite Trapezoidal rule is given without proof in 
Theorem 4.1 


Theorem 4.1 Trapezoidal composite rule 

2 

if f [a, b] , there exists a // e [ a,b ] for which 

Trapezoidal composite rule over n subintervals of 
[( 2 , b\ can be expressed with the error term as 

Pf(x)dx = hf(a)+ 

a 2 

2 znxj)+fm 




Where a=X() < X\ <X2 . < Xyi -1 K Xn = b > 


h=(b-a)/n, and 
For each j=0,l,2,....,n 

Example 2.7 Here , We write a Matlab program to 


approximate fl.6 e X c \ x 
J 1 . 1 


using the Trapezoidal composite rule given by Eqn.(4.1).The 
program is in Fig(4.1) 

gives the exact value , the approximated value and the 
approximated Error by printing values of Ie, Ir and Error 
respectively 


t f(x)dx = ^[f(a) + 2"l f(X2j) 

3 j = l 

m 

+4 Z /QC2/-i) + /(^)] 
7=1 


Ja-b)h 4 ( 4 ) 

180 

Where a= *0 < X\ <X2 . < X2m-\ < X2m 

=b , h=(b-a)/2m, and j=0,l,2,3,....,2m 

Example 2.8 

We compute dx using the 

Simpson’s composite rule given by Eqn.(4.2).The program 
in Fig(4.2) gives the exact value , the approximated value 
and the approximated Error by printing the values of Ie, Ir 
and Error respectively. Observe that Simpson’s rule in 
Example 2.8 uses seven points because We must have 
the number of points to be even and one must be cautious 
about the factors 2&4 in theEqn.(4.2). This why some people 
prefer to apply the trapezoidal rule in Eqn.(4.1) instead of 
Simpson’s rule. 


C. Gaussian quadrature 

All the Newton-Cotes formulas or the formulas that we have 
been given so farrequire that the values of the function whose 
integral is to be approximated be known at evenly spaced 
points, which might be the expected situation if tabulated data 
for the function was being used. If the function is given 
explicitly, however, the points of evaluating the function 
could be chosen in another manner, which leads to increase 
the accuracy of approximation. 

Gaussian quadrature it chooses the points values for 

X\ ,X2 . >X n -1 X n 

in the interval [a ,b ]and the constants 

C\ ’^2 .’ Cyi— l 5 Cfi 

in an optimal manner to minimize the error obtained in 
performing the approximation given by Eqn.(4.3) 

\f(x)dx*> 'Zcjf(xj) ( 4 - 3 ) 

a 7=1 


B. Simpson’s Rule 

Theorem 4.2 Simpson’s composite rule If 

feC\a,b] , there exists a fi e [a, b ] for 

which Simpson’s composite rule over n=2m 
subintervals of \(J,b\ can be expressed with the 
error term as 


Using the Legendre polynomials and their corresponding 
roots , it can be shown quadrature rule is exact for any 
polynomial P(x) ( and has degree of precision at least 
2 n-l)given by 
Eqn.(4.4) . 

1 n 

J p(x)dx — LcjPW (4-4) 
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The roots of Legendre polynomials are the abscissa points 

X\ , X2 > ••••> X n -1 * Xn and the 

coefficients C\,C2 > •• • Cfl —1> both 


are given in Table 4.1 below 

Table 4.1 lists the values for these coefficients and 
abscissa points for n=2,3,4. Other can be found in 
Stroud and Secrest[18] or elsewhere . 

Since the simple linear transformation 
t=[l/(b-a)](2x-a-b) will translate any interval [a, b] 
into [-1, 1] provided b>a , the Legendre polynomials 
Eqn.(4.4) can be used to approximate 


n 


f 


J / {x)dx = Z C j f 
7=1 


a 


(b - a)tb + a 


\ 


v 


(b - a) 


dt 


J 


(4.5) 


for any function that can be evaluated at the required points 


Example 2.9 Approximate the integral 

using Gaussian Quatrature 

rules with n=2 and n=3 .The exact value to seven decimal 

places is 0.1093643 


Table 4.1 shows the coefficients and points for Gaussian rules 

for n=2,3&4 


n 

Roots 

Coefficients 

2 

0.5773502692 

1.0000000000 


-0.5773502692 

1.0000000000 

3 

0.7745966692 

0.5555555556 

0.0000000000 

0.SSSSSSSS89 


-0.7745966692 

0.5555555556 

4 

0.S611363116 

0.347S54S541 

0.3399S10436 

0.6521451549 


- 0.3399S10436 

0.6521451549 


0.S611363116 

0.347S54S541 


% Example 2.7 Integrate using the 
% Composite Trapezoidal Rule 
% In Matlah command window 
% a=l.l; 

% b=l.6; 

% n =5 ; 

% syms x; 

% f=inline ( r exp (x) ' ] ; 

% Ie=compTrapezoidal(f r a,b r n] 


function [Ie P Ir , Error]=compTrapezoidal{f , a,b , n) 
h={b-a]/n. - _ 


x=a:h:b; 
sum=0; 
for i=2:n 

sum=sumffeval { f r x { i ]) ; 
end 


» a=l.l; 
» b=1.6; 
» n=5; 
» s^ttis x; 


% Ie the approximate Integral 
Ie={h/2) * (feval (f P a]+feva.l{f P b)+2*sum) 
% Ir the exact value of Integral 
Ir=int ( r exp{x] r r a,b); 

% The absolute Actual 
% Error = abs{Appoximation - Exact] 
Error = abs(Ie-Ir) r - 


» f=mline('exp(x) ! ) ; 

» 

I e=c ompTr ap ezoi dal(f. ih .n) 
Ie = 1.9505 
Ir = exp(8/5)-exp(l 1/10) 


Fig.(4.1) The Composite Trapezoidal Rule in file 
compTrapezoidal.m & its command window for running it. 


% 4.3 Integrate using the Composite Simpson Rule 
% In Matlah command window 
% a=l . 1 ; 

% b=l . 7 ; 

% n = 3 ; % n is even i.e. n=2*m divisions in the formula 
% syms x ; 

% f=inline { r exp(x) r ] * 

% [Ie,Ir f Error]=compSimpson(f r a r b r n] 


function [Ie,Ir,Error]=compSimpson{f,a,h l n] 
h={b-a]/{2*n] ; 
x=a:h:b r 1 

sum=feval{f , a]; 
for i=2:2:2*n 


Fig,(4.2) (a) The Composite Simpson Rule in 
fil e comp Simps on, m 


sum=sum+4*feval(f,x(i))+2*feval{f r x{i+l]]; 
end 

smn=sum-feval{f ,b); 

% Ie he approximate Integral Ie 
Ie = (h/3) *sran; 

% Ir the exact value of Integral 
Ir=int ( r exp{x] r ,a,b); 

% The absolute Actual Error = abs(Appoximation - Exact] 
Error = abs(le-lrj; 


% 4.3 Integrate using the Commposite Simpson Rule 
% In Matlab command window 
% a=l . 1; 

% b=l . 7; 

% n = 3 ;■ % n is even i.e. n=2*m divisions in the formula 
% syms x; 

% f=inline{ r exp(x] r ] ; 

% [Ie f Ir r Error]=compSimpson(f P a f b,n) 


function [Ie , Ir,Error]=compSimpson(f r a P b,n] 
h={b-a] / (2*n] ; 
x= a :h: b; 

sum=feval{f r a]; 
for i=2:2:2*n 


Fig,(4.2) (a) The Composite Simpson Rule in 
fil e c omp Simp s on.m 


sum=sum+4*feva,l (f r x(i] ] +2*feval (f r x{i+1] ] ' 


end 


sum= sum-f eval { f , b] ; 

% le he approximate Integral Ie 
Ie = (h/3]*sum; 

% Ir the exact value of Integral 
Ir=int ( r exp(x) r ,a P b]; 

% The absolute Actual Error = abs (Appoximation - Exact] 
Error = abs(Ie-Tr); 
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» syms x; 

% Example 4.3(a) Integration using 

» f==inline( l exp(-x 2) 1 ); 

% Guassian Quatrature rules 

» c=[l.0000000000 1.0000000000] 

% In Matlab command window 

» x=[0.5773502692 -0.5773502692]; 

% syms x; 

» a=l; b=1.5 : 

% f=inline(' exp (-x ' 2)'); 

» n=2 ; 

% c=[l.0000000000 1.0000000000]; 

» Ie=quassian(f.c.x.a,b,n) 

% x=[0.5773502692 -0.5773502692];. 

Ie= 0.1094 

% a=l; b=1.5 ; 


% n=2 


% Ie=quas sianif.c.x. a.b.n) 


functi on I e=quas sianf f . c .x. a.b .n) 


sum=0; 


for j=l:n 


t=((b -a)*x(j )+a+b)/2;. 


smn=smn+c(j )* feval(f 3 t)* (b -a)/2 : 

Fig.(4.3(a}} Integration using 

end 

Guassian Quatrature rules in the 

Ie=sum; 

file Quassian.m 


A. Gismalla’s Rules 

Here , We are cited two rules from Gismalla[ 7] , where the 
first one is a fourth degree using seven points to integrate an 
integrand f(x,y) with weight function \'x on the rectangle 
region of integration 

H { (x, _yj : 0 < x < 1 and — 1 < y < 1} . 

The Rule is given by Eqn.(4.6 ), where the 
abscissas x(j) & y(j) and the coefficients c(j) for j=l(l)7 are 
given inTable 4.2 . Similarly the second is a fifth degree rule 
using seven points . The Rule is by Eqn.(4.6), where the 
abscissas x(j) & y(j) and the coefficients c(j) for j=1(1)7 are 
given in Table 4.3 .The Matlab Program for Gismalla fourth 
degree Rule is given in by the file GismallaRule4.m in 
Fig.(4.4) & for the fifth drgree Rule is by GismallaRule5.m in 
Fig.(4.5) 


» syms x; 

» f=inlin ef exp(-x 2 ) 1 ) ; 

» a=l; b=1.5: 

» n= 2 : 

» [I e : Ir : Error]=qua.s sianTable(f .ib ,n) 

Ie = 0.1094 

Ir = 1 fl * erf(3/2)*pi A { 1 12)- 1 !2 * erf (1 )*pi A ( 1 fl) 
Error - 3941559804514209/36028797018963968 
-1 12 * erf(3/2)* pi A ( 1 /2>f 1 fl * erf (1 )*pi A ( 1 12 ) 


Fig.(4.3(b)) The command window for the file 
quassianTable.m with the Gauss Table not passed as 
parameters as in Fig.(4.3(b)) below 


% Example 4.3(b)Integration using Guassian Quatrature rules 

% In Matlab command window 



% syms x; 

% f=inline( ' exp( 
% a=l; b=l.5 ; 

-x A 2)'); 



% n=2 ; % Guassian-Table is given and n runs 

from 2 to 5 

% [Ie , Ir , Error]= 

quassianTable(f,a,b,n) 


function [Ie.Ir.Error] 

=quassianT able(f, 

a.b.n) 


c=[l.0000000000 

0.5555555556 

0.3478548451 

0.2369268850 

1.0000000000 

0.8888888889 

0.6521451549 

0.4786286705 

O.0000000000 

0.5555555556 

0.6521451549 

0.5688888889 

0.0000000000 

0.0000000000 

0.3478548451 

0.4786286705 

O.0000000000 

0.0000000000 

0.OOOOOOOOOO 

0.2369268850]; 

x=[ 0.5773502692 

0.7745966692 

0.8611363116 

0.9061798459 

-0.5773502692 

0.0000000000 

0.3399810436 

0.5384693101 

O.OOOOOOOOOO 

-0.7745966692 

-0.3399810436 

0.OOOOOOOOOO 

O.OOOOOOOOOO 

0.0000000000 

-0.8611363116 

-0.5384693101 

0.0000000000 

0.0000000000 

0.0000000000 

-0.9061798459]; 

sum=0; 
for j=l:n 




t=((b-a)*x<j,n-l)+a+b)/2; 
sum=sum+c(j,n-l)*feval(f,t)*(b-a)/2; 



end 

Ie=sum; 




Ir=int('exp(-x A 2)'.a,b); 



Error = abs(Ie-Ir); 





Fig.(4.3(b)) Guassian Quatrature rules in the file 

quassianT able.m 


V . Cubature ’ s Rules for 2-dimentional Integration 

Cubature formulae are quadrature formulae that having 
the minimum points to attain the required degree of accuracy 
when evaluating a polynomial as an integrand .In a series of 
papers rules for constructing cubature formulae have been 
established as in STRUOD [4], SCHIMD [17],RADON[14] 
andGISMALLA[ 6 - 8 ]. 


i i 



: -l 


f(x,y) dxdy ® 


j=i 



B. Radon’s Rule 


Radon[ 14 ] has established a fifth degree formula where the 
region of the integration is the square with the unity weight 
function w(x)=l and the formula is given by Eqn.( 4.7) 


Fig(4,4) The fou. ith degret Gh m allaRulel with 

Madab Program 

GhmaUaRuMm & irs Comm a ad Window 


% Integration using GismailaRule for Double Integration on a rectangle . 
% This 'Quatrature rule of fourth degree to evaluate an intgrand 
% x A 0.5*f(s,y) on the regin R= {{x : y}: 0=< x =<1 and - 1 =< y =< 1 } 

% The integrand must be called without the weight function sqrt(x) 

% as an inline or function handle object of the form f(x ; y) only. 

% In Matlab command window 
% syms s; 

% syms t; 

% f=mline( B exp{-s*t) r ); 

% a= 0 ;b=l: 

% c=-l;d=l ; 

% n=7; % Gismalla-Table is given for n = l 
% [Ie]=GisniallaMe4(f,a,b, c,d,n) 


function [Ie]= Gi smallaRule4{f. a. b. c.d. n) 
x=[ 0.0607124836 0.0000000000 0 . 
0.3323846207 0.5329336163 0. 
0 1 .19949548 83 
0.1851851852 
0.1851851852 
0.1851851852 
0.1851851852 
TotalSum=0; 
for j=l:n 

5(j)={{b-a) , "x{j ! 2)+a ) ; 
t(j)=((d-c) *x(j,3)+c+d)/2; 
sum{j)= x(j, l}*feval(f, s(j) , t(j) ) ; 
TotalSmn=TotalSum - sum{j); 

End 

I«={d-c)*0.5*TotalSum : 



0.8943391109 0. 

0.8618614683 0.7745966692 
0.8618614683 -.7745966692 
0.3381385317 0.7745966692 
0.3381385317 -.7745966692 ]; 


» ?. : ms s; 
s™ t; 

j * 

»f=mlm<' e?:p(-ih) r ); 

ap=0 ; t=l ; 

» c=-l; d=l ; 

» n=7; 

» |Te]= GismallaRule4(f.a.b s c. d .n) 1 
» fc=1.4J18 
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Table 4.2 The abscissas x(j) & y(j) and the coefficients c(j) 
_ for j=l(l)7 for the fourth degree _ 


j 

c(i) 


vU> 

I 

0.0607124336 

#..########## 

0,0000000000 

2 

0,3323346207 

0.53293345163 

0.0000000000 

3 

0,1994954333 

0,3943391109 

0,0000000000 

4 

0,1351351352 

0,3613614633 

0.7745966692 

5 

<0,1351351352 

0,3613614633 

0.7745966692 

6 

O. 1351351352 

0.3331335317 

0.7745966692 

T 

#.1351351352 

0.3331335317 

0.7745966692 


Table 4.3 The abscissas x(j) & y(j) and the coefficients c(j) 
for j=1(1)7 for the fifth degree 



cC) 


Vft> 

J 

9.3911751535 

9.5£S£1$5251 

9. 9999999999 

2 

9.1395S37632 

9.9$719$£2££ 

9. 9999999999 

3 

9. 979$33£75£ 

9. 9596574637 

9. 9999999999 

4 

9.154197776$ 

9.2$99491979 

9.77459£££?2 

5 

9.154197776$ 

9.2$99491979 

9.77459£££?2 

£ 

9.21 £1725936 

9.$211£19132 

9.7745966692 

spy 

-s 

9.21 £172593£ 

9. $211619132 

9.7745966692 


% Integration using Gismalia Rule for Double Integration on a rectangle . 
% This Quatrature rule of fifth degree to evaluate an intgrand 
% x J,,, 0.5*f(x,y) on the regin R={(s 3 y}: 0=< x =<1 and -1=< y =<1} 

% The integrand must be called without the weight function sqrt(x) 

% as an inline or function handle object of the form f{x.y) only. 

% In Matlab command window 
% symss; 

% svmst: 

% f=inline{" exp(-s*t)); 

% a=0; b=l ; 

% c=-l; d=l : 

% n=7; % Gismalia-Table is given for n = 7 
% [Ie]= GismaliaRule5{f : a : b F c, d m) 


function [Ie]= Gi smallaRule5{f. a.b.c.d. n) 
x=[ 0.391175153S 0.5686185251 0.0000000000 
0.1305837632 0.9871086266 0.0000000000 
0.0708336756 0.0596574637 0.0000000000 
0.1541977768 0.2899491979 0.7745966692 


0.1541977768 0.2899491979 
0.2161725936 0.8211619132 
0.2161725936 0.8211619132 
Total Sum=0; 
forj=l:n 

^)=0>a)*x(j,2)+a ) ; 
t(j)K(d-c) *x(j,3)+c+d)/2; 
sum{j)= x(j, l)*feval(f a sQ), t(j) ) 
TotalSmn=Total Sum h 

End 

Ie={d-c)*0.5*TotalSum : 


-.7745966692 

0.7745966692 

-.7745966692 


]; 


» syms s: 

» syms t; 

» f= tnltn rt T COi(0.:^ : pb : (3+t)) r >; 

» a=0; b=l ; 

» c=-l ; d=l 
---■ u=7; 

» [Ie]= Gi3ni=U=RuIe5{f a,b. c, d .n) 
» Ie= 1.4316 


Fig(4.5) The fifth degree GismallaRule5 with Matlab 

Program 

GismallaRule5.m & its Command Window 


l l 


| | fCx, y) dxdy w ^ c(j3f(xCj3,yCp) 

-i -i i=i 



Where the abscissa and weight coefficient are in Table 4.4 
and the Matlab Program for RadonRule.m is in the Fig.(4.6). 
As in Gismalla[ 12 ] , the reader can test these the three 
programs GismallaRule4.m, GismallaRule5.m and Radon 
Rule 5.m for the three given integrals below with their exact 
values 


i i 


-// 


1 1 


Vx cos [ — (x + y) | dxdy ^ 0.45 5 33 (4. B) 


o -1 


i l 



1 1 


and ], = 


L = I I Vx e-^'dxdy w 1.42196496 (4.9) 

Z -1 

1 


£ 

VI 


J J 


mv :v 


I! -1 


4+jf+y 


0.1956196 ( 4,10 


% Integration using RadonRuIe for Double Integration on a square . 

% This Quatrature rule of fifth degree to evaluate an intgrand 
% f(x : v) on the re gin R={(x.y): -1=< x =<1 and -1 =< y =<1} 

% The integrand must be called with the unity weight function 
% as an inline or function handle object of the fomif(x : y) . 

% In Matlab command wmdow 

% symss; Tike fifth Rfl-dauRulet niflk Matlab Progrui 

Vo svmst; 

% f=inline( sqrt(s)* exp(-s + t) r ); 

%a=G;b=l; 

% c=-l ; d=l ; 

% n= 7 ; % Ra don-Table is giv en for n = " 

% [Ie]= RadonRuIe 5(f a.b. c.. d,n) 

function [Ie]= RadonRuIe5(f.a,b,c,d.n) 
x=[ S/7 0 0 

20/63 0 (14/15/0.5 

20 63 0 -<14/15/0.5 

5/9 (3/5/0.5 (1/3/0.5 

5/9 -(3/5/0.5 (1/3/0.5 

5/9 (3/5/0.5 (1/3/0.5 

5/9 (3/5/<)_5 4:l/3/0.5 ]; 

TotalSum=0; 
forj=l:n 

s(iH(b-a)*x(j,2)fa+by2 ; 
t(j ;;K(d-c)*x(j ,3 )+c+d yi ; 
sum(j )=x(j4 )* fe v alffsfj )4(j)); 

T ot ,alSum=T ot alSum- sum(j); 
end 

TotalSum= (b-a)*(d-c)*TotalSum 4; 

Ie=TotalSum ; 


» syms s; 

» svmst; 

» f=inline( T £qrt(s)*exp(-s H: t) f ); 

» a=0;b=l ; 

» c=-l; d=l ; 

» n=7; 

» Pe]= RadonRuIe5(f.,a,b. c.. d ,n) 
» Ie=l .43 563 SS622271S 


Table 4.4 The abscissas x(j) & y(j) and the coefficients c(j) 
for j=l (1)7 for RadonRule5 in Eqn.(4.7) with Matlab 


Program RadonRule5.m in Fig(4.6) 


■ 

I 

*U) 

iG) 

>0) 

I 

8/7 

0 

0 


20.63 

0 

(14/15) "0.5 

3 

20/63 

0 

- (14/15) M). 5 

4 

5/9 

(3/5) M). 5 

(1/3) J -0 .5 

8 

5/9 

-(3/5) M).. 5 

(1/3) *0.5 

6 

5/9 

( 3/5) '-0.5 

(1/3) '-0 .5 

7 

5/9 

-(3/5) .5 

- ( 1/3) MS .5 


VI. Levin’s Transform Rule for higher dimensional 

Integration 

Levin’s U-transform Method [12] has been shown and 
used to efficiently compute certain types of multiple integrals 
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to a high decimal places of accuracy using FORTARN 
program languages . The Computing Machine was the main 
FRAME computer applied with a FULL machine accuracy to 
get high decimal places of accuracy. The Levin’s Method has 
a drawback; it fails to compute certain types of series. 

Instead of FORTARN program language written for 
Levin’s Method which cannot be available now or after, We 
write it here instead using Matlab program. In contrast for 
using FULL machine accuracy when using FORT AN 
language, the command format long can be used instead in 
the Matlab program to get high decimal places of accuracy . 
Levin’s U-transform is defined in Eqn.(4.11) as 
U 2 #k+i_(X) = 

s j + i 


Zk+l- 


Ej = T(-i) ] c J 


)G+i) zk 

+ 1 




(+. 11 ) 


where S r . is the partial sum of n terms of the convergent series 
of positive terms as in Eqn.(4.12). 



l 


Here , our aim is to apply Levin’s U-transform Matlab 
program to evaluate the three multiple integrals considered in 
[12] using FORT ARAN language . The three Lattice Green 
functions are 

1 r T 

G(0Q0) = — f 

7T J jj, 


rr _ t 

J c 1 - CO. 


ax ay az 


s x cos v cos z 


(4.13) 


G( 100 ) = 




cos 2 x ax ay c'z 
1 — cos xcos y cos z 


(4.140 



cos2xcos2y dxdy dz 
1 - cosjccosy cosz 




Since cos(2x)=cos 2 (x)-l and 

cos(2y)=cos 2 (y)-l , it follows that 
GtlOti) =2I 4 - GiOOO 0 (4.16) 

GO 10) = 4I 3 - 41 4 + G{QQQ) (4.17) 


where 


?t -7T pn cds*x dx dy ds 




C D 8 X C D £y C D 


(4.IS) 


l r r 1 f cos 1 . 

Ij ~ ^ J E J J "TT 


2 x cos 2 x dx dy dz 


cos x cos y cos z 


Now by expanding the integrand in G(000), 1 4 and 1$ in 
Eqn.(4.13) ,Eqn.(4.18) and Eqn.(1.19) respectively and 
applying Wallis’ formula , We get the integrals as 

■E 

G(000) = Ya£ (4.20) 

<n 

'4 = 2 ^^+! C 4 ’ 21 ) 

I 5 = sEo a r.fl 0.22) 


Table 4.5 LEVIN TKAN5F OR} 1 SUM FOR G{000) 


k 

The- Sum 

LE\T> TRAN5F OR} lU 2mk ^ 

0 

1.00000000 

1.4000720 

I 

1.12500000 

1.3931052 

■j 

1.17773437 

1.39320476 

3 

1.20825195 

1.3932039225 

4 

1.2286904 

1.39320392963 

5 

1.24360030 


5 

1.25508015 


7 

1.26427156 


B 

1.27184504 


9 

1.27822511 

G{000>= 1.3932039296B 

10 

1.28369522 


11 

1.28845279 



Table 4 6 LE\T>~ TRANSFORM 1 SUM FOR l 4 


k 

The-Sum. S k LEVIN TKANS’OKU + 

0 

0.500000000 

0.740888952 

1 

0,593 750000 

0,832874813 


0.637695312 

0.841755935 

3 

0, 66439B193 

0,842047393 

4 

0,6 B2798147 

0,842052522 


0,696460112 

0,842052578 

6 

0.70 7119970 


7 

0,715736914 


8 

0.722889651 

l 4 =0.842052578 

9 

0.728950713 


10 

0.734172180 

G(1O0)=0 290901226 

11 

0,738731524 


12 

0.742757787 


13 

0.746347348 


Table 4.7 LEATN TKANSFOKA1 SUM FOR L 

k 

The- Sum S k LEVIN TELANSF OKA I U 2mk4 . : 

0 

0.250000000 

0.380674380 

1 

0.320312500 

0.533267522 

■J 

0.356933593 

0.550563530 

3 

0.380298614 

0.551141114 

4 

0.396B5B572 

0.551151238 

b 

0.409382041 

0.551151349 

6 

0.4192804 BG 


7 

0.427358865 


8 

0434114228 1 

[-= 0551151349 

9 

0.439872237 


10 

0.444856364 

G(100) =0.229599014 

11 

0.449225735 


12 

0.453097143 


13 

0.4565 58504 
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% Integration using Levin Transform for Triple Integration on cubic. 

% This Technique is a series Technique for which the integral is expressed 
% as a series FIRST and then LEVIN is applied. The general term for the 
% series can be expressed as an inline function or a handle object with 
% the initial term submitted to the program in advanced to generate the 
% other terms with the number of terms n to be taken for the sum. 

% The function fin LevinTranfSum(f : n : a0) is a ratio to generate other terms . 

% In Matlab command window 
% syms s; 

%f=inlineC ((2*s-l)/(2*s)) A 3'): 

% aO=l ; 

%n=5; 

% LevinTranfSum(f ; n : a0) 
function LevinTranfSum(f : n : a0) 
global UT; 
for k=l:n 

[S : UT]=Le\inTransform(Tk : a0); 

F(k)=UT; 
end 

dispC The Sum of the Series having 2*n+2 terms') 
disp([ S']): 

disp(' The Sum of the Series using LEVIN TRANSFORM using 2* n+2 terms') 
disp([ F’]): 

function [ S , UT]=LevinTransform(f : k : a0) 
a(l)=aO ; 


» format long 
»svms s: 

» f=inlineC ((2*s-l)/(2*s)) A 3'); 
»aO=l : 

»n=5; 

» LevinTranfSum(f : n : a0) 


The Command Window- for G(000) and 
its results are in Table (4.5) 


S(l)= a(l); 

C(l)= 1; 

TotalSumDen(l)=l; 
TotalSumNum(l)=l: 


for j=l:2*k+l 

a(j+1 )=feval(fj)*a(j); 

S(j+l)=S(j)+a(j+l); 

C(j+lH2*k+2-j)*C(j)'(j); 


> > format long 
»syms s: 

» f=inlineC ((2*s-l>'(2*s)) A 2*((2*s+l)/(2*s+2))'); 
» a0=0.5 ; 

» n=6; 

» LevinTranfSum(f : n : a0) 

The Command Window- for I 4 and 
its results are in Table (4.6) 


TotalSumDen(j+l)=TotalSumDen0 + (-l) A j*C(j+l)*(j+l) A (2*k-l)*S(]+l)/a(j+l); 
TotalSumNum(j+l)=TotalSumNum(j) + (-l>"j*C(j+l)*(j+l) A (2*k-l), , a(j+l); 


end 


UT= Tot al SumD en(2 *k+2) T ot al SumN tun (2 *k+2); 


Fig(4.7) Matlab Program LevinTranfSum.m with its Two 

Command Window 

for G(000) and [4 . The Output are in Table 4.5 &Table 4.6 


where the coefficients terms a’s are given by 


1.3.5 . (2n - 1) 

2..4.6 ... (2nj 





Now the Matlab Program to integrate G(000) , 1 4 and in 
Eqn.(4.20) ,Eqn.(4.21) and Eqn.(4.22) respectively is given in 
Fig(4.7). The output 

results from the command window program for G(000) , 1 4 
and /j are given in Table 4.5 , Table 4.6 and Table 4.7 
respectively.. 

Observe that to compute the integral I 5 , 

We need to replace the generating ratio function 

and the initial term aO by 

f=inline(' (( 2 *s-l)/( 2 *s))*(( 2 *s+l)/( 2 *s+ 2 )) A 2 '); 

and a0=0.25 ; in the command window program for 

LevinTranfSum.m given in Fig.(4.7) . 

Hence the values for the integral G(100) in Eqn.(4.16) 
and G(110) in Eqn.(4.17) can be 
evaluated up to eight decimal places of 
accuracy as in Eqn.(4.24) where 
G(110)=0.229599014 

and G( 100)=0.29091226 (4.24) 


VIE Conclusion 

The reader must know all types of Programs Built-in as 
ROUTINES SOFTWARE 


are inflexible programs to suit and solve any particular types 
of problems while the others are flexible programs. This can 
be seen if someone attempts to run the routine DBLQUAD on 
MALAB for the integral 

I ] in Eqn. (4 . 1S} , U m Eqn. (4.1 9 ) 

or I 3 in Eqn. (4.20 j an ERROR will be occurred. This 
because DBLQUAD on MATLAB for interval [0,1] which 
runs on variable x will be transformed for a different one by 
the inside QUAD quadrature used by DBLQUAD . Even if 
one doesn’t submit the parameters a =0 and b=l in the 
MATLAB program GismallRule4.m or GismaleRule5.m an 
ERROR will be occurred due to the fact that any 
transformation to the interval [ 0 , 1 ] will effect to translate the 
corresponding weight function Vx given in Eqn, (4.6 j to a 
different weight function and so the abscissas’ and 
coefficients will be different from their corresponding in the 
given rule in Eqn.(4. 6 ) 

However, the program RadonRule5.m can run these 
integrals efficiently without any errors provided the integrand 
function must be submitted as 

\ x fix, yj exactly as with these 
integrands given by Eqn.(4.18), Eqn.(19) and Eqn.(20). 

Levin’s Transform in Eqn.(4.11) written as functions 
or subroutines using FORTRAN languages , even on 
main-frames computers will compute the integrals G(000), 1 4 
and/jin Eqn.(4.20) ,Eqn.(4.21) and Eqn.(4.22) respectively 
only to five or sixth decimal places at most. In Gismalla[12] 
these integrals were computed to a very high decimal places 
up to 12 decimals places using the same Levin’s Transform in 
Eqn.(4.11) . The cost is very very too expensive , the 

technique Levin’s Transform must be phrased and 
unexpressed as function or procedures each times when it is 
called with SOME SPECIAL COMMAND WITTEN ON 
THE FIRST TOP LINES to generate the FULL MACHINE 
ACCRACCY which are unavailable unless given by the 
advisory of the main-frame 
computers. The computational remarks with 
Levin’s Transform as a function written in FORTRAN 
language can found in gismalla[ 12 ]. 

On the contrary, the LevinTrasfSum.m given in Fig.(4.7) 
computes these integrals neatly 

up to eight digits of accuracy exactly without a great 
complexity and efforts using the symbolic languages 
MATLAB and even the program is SHORT. 
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